function [A,xmi,xma,inx,ymi,yma,iny,nr,nc,filenout] = gs_dtmread(filena)
%
%GS_DTMREAD Read a z-matrix related to a DTM generated by Surfer (ASCII 
% .grd file) or ArcGIS (ASCII .asc file). Syntax:
%
%[A,xmi,xma,inx,ymi,yma,iny,nr,nc,FILENOUT] = GS_DTMREAD,
%
% where A is the output MATLAB matrix (z values assigned on a regular xy-grid).
% The ASCII filename is interactively asked.
% 
% Surfer raster case:
% --------------------
%
% An ASCII grid file by Surfer (.grd) contains five header lines that 
% provide information about the size and limits of the grid, followed by 
% the list of Z values. The fields within ASCII grid files must be space 
% delimited. The general format of such an ASCII grid file is:
%
% ID
%   The identification string DSAA that identifies the file as 
%   an ASCII grid file.
%
% nc nr
%   nc is the integer number of grid lines along the X axis (columns)
%   nr is the integer number of grid lines along the Y axis (rows)
% 
% xmi xma
%   xmi is the minimum X value of the grid
%   xma is the maximum X value of the grid
%
% ymi yma
%   ymi is the minimum Y value of the grid
%   yma is the maximum Y value of the grid
%
% zmi zma 
%   zmi is the minimum Z value of the grid
%   zma is the maximum Z value of the grid
%
% grid row 1
%   ...         these are the rows of Z values of the grid, organized in 
% grid row nr   row order. Each row has a constant Y coordinate. So, grid
%               row 1 corresponds to ymi and the last grid row corresponds
%               to yma. Within each row, the Z values are arranged from 
%               xmi to xma. This agrees withe the matlab representation
%               (in particular, FLIPUD function is not necessay)
%
% ESRI ArcGIS raster case:
%--------------------------
% 
% An ESRI ASCII grid file by ArcGIS [.ASC] contains six header lines that 
% provide information about the size and limits of the grid, followed by 
% the list of Z values. The fields within ASCII grid files must be space 
% delimited. The header of an ASCII .asc grid file is:
%
%   ncols           xxx
%   nrows           xxx
%   xllcorner       xxx
%   yllcorner       xxx
%   cellsize        xxx
%   NODATA_value    xxx
%
% or, alternatively,
%
%   ncols           xxx
%   nrows           xxx
%   xllcenter       xxx
%   yllcenter       xxx
%   cellsize        xxx
%   NODATA_value    xxx
%
% where
%   
%   ncols (respectively nrows) is the integer number of grid lines along 
%   the X-axis (Y-axis)
% 
%   xllcorner (yllcorner) is the X-coordinate (Y-coordinate) of the grid 
%   origin by lower-left corner of the cell.
%   Alternatively, xllcenter (yllcenter) is the X-coordinate (Y-coordinate)
%   of the grid origin by center of the cell.
%
%   cellsize is the cell size (rqual for X- and Y- direction)
%
%   NODATA_value is the value related to no data cells (NaN in the oputput
%   matrix). Although the row is compulsory (the header or an ESRI raster
%   file must have six rows), the value can be omitted. In this case, the
%   default NODATA_value is considered, i.e. -9999.
%   
% grid row 1
%   ...         these are the rows of Z values of the grid, organized in 
% grid row nr   row order, starting from the bottom, left value. Therefore,
%               FLIPUD function is automatically called in order to
%               correctly align the data coording to MATLAB conventions.
%
% Other output variables:
%
% The output cells sizes are inx and iny respectively (inx and iny can be
% different in the case of a SURFER-like file, but have the same value for
% and ESRI file.
% The string FILENOUT is the name of the opened .grd or .asc file. 
%
%[A,xmi,xma,inx,ymi,yma,iny,nr,nc,filenout] = GS_DTMREAD(FILENA).
%
% The function operates as above, but the data are read from the file whose
% filename is FILENA (.GRD or .ASC file, case insentitive for the extension).
% Also in this in this case, FILENOUT is the effectively opened .grd or 
% .asc file (if FILENA is opened, it is FILENOUT = FILENA). If FILENA is 
% invalid, the filename is interactively asked).
%
%   [A,xmi,xma,inx,ymi,yma,iny,nr,nc,filenout] = GS_DTMREAD
%   [A,xmi,xma,inx,ymi,yma,iny,nr,nc,filenout] = GS_DTMREAD(FILENA)
%
% See also GS_DTMWRITE.

% G. Teza, 2006, 2015.

if nargin == 0, filena = []; end

nowell = 0;
k = 0; 
while nowell == 0
    
    if isempty(filena)||(k > 0)
        [filen, pathn] = uigetfile(...
            {'*.grd','Surfer ASCII files (*.grd)';...
            '*.asc','ESRI ArcGIS ASCII files (*.asc)'},...
            'SURFER (.GRD) OR ESRI (.ASC) DTM ASCII FILE');
        if filen
            filena = fullfile(pathn,filen); 
        else
            filena = [];
        end
    end
    
    if ~isempty(filena) && strcmpi(filena(end-2:end),'grd')
        
        try 
            nind = dlmread(filena,' ',[1 0 4 1]);  % extraction of parametrs
            nind = reshape(nind',1,8);       % a row-vector must be the result! 
            nc  = nind(1); nr  = nind(2);
            xmi = nind(3); xma = nind(4);
            ymi = nind(5); yma = nind(6);
            zmi = nind(7); zma = nind(8);
            inx = (xma-xmi)/nc;
            iny = (yma-ymi)/nr;
            if abs(inx-iny) < 1000*eps, iny = inx; end
            idrow1 = 5; % index of initial row for DLMREAD function
            nowell = 1;
        catch    %#ok<CTCH>
            if ((nargin == 1)&&(k > 0))||(nargin == 0)
                menwar = menu(...
                    sprintf('THE .grd ASCII FILENAME %s IS INVALID',filena),...
                    'ANOTHER FILENAME','EXIT');
                if menwar == 2, nowell = 1; filena = []; end
            end
        end
        
    elseif ~isempty(filena) && strcmpi(filena(end-2:end),'asc')
        
        try
            fid = fopen(filena);
            C = textscan(fid,'%s%f');
            fclose(fid);

            SC  = C(1);  % string values
            SCD = deal(SC{:});
            NC  = C(2);  % numerical values
            NCD = deal(NC{:});

            nc = NCD(1);  nr  = NCD(2);
            inx = NCD(5); iny = inx;
            if strcmpi(char(SCD(3)),'xllcorner')
                xmi = NCD(3);
            else
                xmi = NCD(3)-dx/2;
            end
            if strcmpi(char(SCD(4)),'yllcorner')
                ymi = NCD(4);
            else
                ymi = NCD(4)-dx/2;
            end

            noval = NCD(6);
            if isnan(noval), noval = -9999; end % in this case, the 6th header value is undefined
            xma = xmi+nc*inx; 
            yma = ymi+nr*iny;
            idrow1 = 6; % index of initial row for DLMREAD function 
            nowell = 1;

        catch %#ok<CTCH>
            if ((nargin == 1)&&(k > 0))||(nargin == 0)
                menwar = menu(...
                    sprintf('THE .asc ASCII FILENAME %s IS INVALID',filena),...
                    'ANOTHER FILENAME','EXIT');
                if menwar == 2, nowell = 1; filena = []; end
            end
        end
            
    else
        menwar = menu(...
            'INVALID FILENAME. IT MUST BE A .grd (Surfer) OR AN .asc (ArcGIS) FILE',...
            'ANOTHER FILENAME','EXIT');
        if menwar == 2, nowell = 1; filena = []; end
    end
    k = k+1;
end
if isempty(filena) 
    A = []; 
    xmi = []; xma = []; inx = 0; ymi = []; yma = []; iny = 0;
    nr = 0; nc = 0; filenout = []; 
    return 
end

A  = dlmread(filena,' ',idrow1,0); % 5 --> 6th row, 0 --> first column .grd (6 --> 7th row, 0 --> firts column .asc) 
[~,ncl] = size(A);

if ncl == nc+1              % in this case, dlmread introduced an extracolum to be removed 
    A(:,end) = [];
end

if strcmpi(filena(end-2:end),'grd')  
    A = gs_subnan(A,zmi,zma);   % The unmodeled points are transformed in NaN
else
    J = A == noval;
    A(J) = NaN;             % The unmodeled points are transformed in NaN
    A = flipud(A);          % Please note that FLIPUD is necessary if a .asc file is used
end

filenout = filena;